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An ensemble of random unistochastic (orthostochastic) matrices is denned by taking squared 
moduli of elements of random unitary (orthogonal) matrices distributed according to the Haar mea- 
sure on U(N) (or O(N), respectively). An ensemble of symmetric unistochastic matrices is obtained 
with use of unitary symmetric matrices pertaining to the circular orthogonal ensemble. We study 
the distribution of complex eigenvalues of bistochastic, unistochastic and orthostochastic matrices 
in the complex plane. We compute averages (entropy, traces) over the ensembles of unistochastic 
matrices and present inequalities concerning the entropies of products of bistochastic matrices. 



I. INTRODUCTION 

Consider a square matrix B of size N containing non-negative entries. It is called stochastic if each row sums to unity 
(J2iLi Bij — 1 for j — 1, . . . , N) (for information about properties of such matrices consult QJ^]). If, additionally, 
each of its columns sums to unity, i.e., Bji — 1 f° r i = I, ■ ■ ■ , N , it is called bistochastic (or doubly stochastic). 

Bistochastic matrices emerge in several physical problems. They are used in the theory of majorization [|3|-|6|, 
angular momentum j7|, in transfer problems, investigations of the Frobenius-Perron operator, and in characterization 
of completely positive maps acting in the space of density matrices ||. We shall denote by fl% (resp. f2jy) the sets 
of stochastic (resp. bistochastic) matrices of size N. 

A matrix B is called orthostochastic, if there exists an orthogonal matrix O, such that Bij = O, 2 for i, j = 1, . . . , N. 
Analogously, a matrix B is called unistochastic (unitary-stochastic)^, if there exists a unitary matrix U, such that 
B^ = \Uij \ 2 ior i,j — 1, . . . , N. Due to unitarity (orthogonality) condition every unistochastic (orthostochastic) matrix 
is bistochastic. These four sets of nonnegative matrices arc related by the following inclusion: £1^ C C fijy c Slff, 
where QjL and £1% represent the sets of unistochastic (orthostochastic) matrices. For TV > 2 all three inclusions are 



proper 

Unistochastic matrices appear in analysis of models describing the time evolution in quantum graphs |9 13 and 



in description of non-unitary transformations of density matrices [ |14|jl5| . Moreover, the theory of majorization and 
unistochastic matrices plays a crucial role in recent research on local manipulations with pure states entanglement 
fl6| or in characterizing the interaction costs of non-local quantum gates [JL7| . 

In this work we analyse the structure of the set of bistochastic (unistochastic) matrices of a fixed size and investigate 
the support of their spectra. Knowledge of any constraints for the localisation of the eigenvalues of such a matrix 
is of a direct physical importance, since the moduli of the largest eigenvalues determine the rate of relaxation to 
the invariant state of the corresponding system. We define the notion of entropy of bistochastic matrices and prove 
certain inequalities comparing the initial and the final entropy of any probability vector subjected to a Markov chain 
described by an arbitrary bistochastic matrix. A related inequality concerns the entropy of the product of two 
bistochastic matrices. 

Moreover, we define physically motivated ensembles of random unistochastic matrices and analyse their properties. 
As usual, the term ensemble denotes a pair: a space and a probability measure defined on it (for example the circular 
unitary ensemble (CUE) represents the group U(N) of unitary matrices of size N with the Haar measure fL8[ ). Since 
any unitary matrix determines a unistochastic matrix, the Haar measure on the unitary group U(N) induces uniquely 
the measure in the space of unistochastic matrices fijy, which we shall denote by Analogously, we shall put \xq 



1 This notation is not unique: in some mathematical papers (e.g. G3) ) unistochastic matrices are called orthostochastic 



1 



for the measure in fiS induced by the Haar measure on the orthogonal group O(N). In the sequel we shall use the 
names unistochastic ensemble (resp. orthostochastic ensemble) for the pair {ri^,/Z[/} (resp. {Q° : fj,o})- 

We compute certain averages with respect to these ensembles. Related results were presented recently by Berkolaiko 
fl9f and Tanner jjjfl . In the latter paper the author defines unitary stochastic ensembles, which have a different 
meaning: they consist of unitary matrices corresponding to a given unistochastic matrix. 

This paper is organised as follows. In section II we review properties of stochastic matrices. In particular we 
analyse the support of spectra of random bistochastic and stochastic matrices in the unit circle. In section III some 
results concerning majorization, ordering, and entropies of bistochastic matrices are presented. In particular we prove 
subadditivity of entropy for bistochastic matrices. Section IV is devoted to the ensembles of orthostochastic and 
unistochastic matrices; we investigate the support of their spectra, compute the entropy averages, the average traces, 
and expectation values of the moduli of subleading eigenvalues. Some open problems are listed in section V. Analysis 
of certain families of unistochastic matrices and calculation of the averages with respect to the unistochastic ensemble 
is relegated to the appendices. 



II. STOCHASTIC AND BISTOCHASTIC MATRICES 



A. General properties 



In this section we provide a short review of properties of stochastic and bistochastic matrices. The set of bistochastic 
matrices of size N can be viewed as a convex polyhedron in R N . There exist Nl permutation matrices of size ^V, 
obtained by interchanging the rows (or columns) of the identity matrix. Due to the Birkhoff theorem, any bistochastic 
matrix can be represented as a linear combination of permutation matrices. In other words the set of bistochastic 
matrices is the convex hull of the set of permutation matrices. By the Caratheodory theorem it ispossible to use 
only TV 2 — 1 permutation matrices to obtain a given bistochastic matrix as their convex combination |3[] . Farahat and 
Mirsky showed that in this combination it is sufficient to use N 2 — 2N + 2 permutation matrices only, but this number 
cannot be reduced any further |f20f . The dimension of the set of bistochastic matrices is (N — l) 2 . The volume of the 
polyhedron of bistochastic matrices was computed by Chan and Robbins ]2l| ] . 

Due to the Frobenius-Perron theorem any stochastic matrix has at least one eigenvalue equal to one, and all others 
located at or inside the unit circle. The eigenvector corresponding to the eigenvalue 1 has all its components real and 
non-negative. For bistochastic matrices the corresponding eigenvector consists of N components equal to 1/7V and 
is called uniform. A stochastic matrix S is called reducible if it is block diagonal, or if there exists a permutation P 
which brings it into a block structure, 



S' = psp- 



A 1 
C A 2 



(2.1) 



where Ai are square matrices of size Ni < N for i = 1, 2, N = N% + N2. It is called decomposable if one can find 
two permutation matrices P and Q such that PSQ has the above form. Matrix S is irreducible (indecomposable) if 
no such matrix P (matrices P and Q) exists (exist) |3|,^2| . An irreducible stochastic matrix cannot have two line arly 
independent vectors with all components nonnegative. Any reducible bistochastic matrix is completely reducible j22|, 



i.e., the matrix C in (2.1) is equal to zero. 

A stochastic matrix S is called primitive if there exists only one eigenvalue with modulus equal to one. If S is 
primitive, then S k is irreducible for all k > 1 ||. Note that the permutation matrices P with trP — are irreducible, 
but non-primitive, since P N equal to identity is reducible. For any primitive stochastic there exist a natural number 
k such that the power S has all entries positive. The fact that all the moduli of eigenvalues but one are smaller than 
unity implies the convergence lim^oo B k — B*. Here £>» denotes the uniform bistochastic matrix (van der Waerden 
matrix) with all elements equal to 1/N. Its spectrum consists of one eigenvalue equal to one and N — 1 others equal 
to zero. The matrix B* saturates the well known van der Waerden inequality ||] concerning the permanent of the 
bistochastic matrices: peri? > N\/N N , and hence is sometimes called the minimal bistochastic matrix [p3| , p2[ . 

Each bistochastic matrix of size N may represent a transfer process at an oriented graph consisting of N nodes. If 
a graph is disjoint or consists of a Hamilton cycle (which represents a permutation of all N elements), the bistochastic 
matrix is not primitive, and the modulus of the subleading eigenvalue (the second largest) is equal to unity. 

If matrices A and B are bistochastic, its product C — AB is also bistochastic. However, the set of bistochastic 
matrices does not form a group, since in general the inverse matrix A -1 is not bistochastic (if it exists). For any 
permutation matrix P its inverse P _1 = P T is bistochastic and the eigenvalues of P and P _1 are equal and belong 
to the unit circle. 
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B. Spectra of stochastic matrices 



A stochastic matrix contains only non-negative entries and due to the Frobenius-Perron theorem its largest eigen- 
value is real. This leading eigenvalue is equal to unity, since its spectral radius is bounded by the largest and the 
smallest sum of its rows, all of which are equal to 1. In the simplest case of permutation matrices the spectrum 
consists of some roots of unity. The eigenvalues of permutation matrices consisting of only one cycle of length fc are 
exactly the fc-th roots of unity. 

Upper bounds for the size of the other eigenvalues are given in |24|] . Let M denote the largest element of a stochastic 
matrix and m the smallest. Then the radius r 2 of a subleading eigenvalue satisfies 



r 2 < 



M 



M 



From this bound it follows that all subleading eigenvalues of the van der Waerden uniform matrix 
Another simple bound of this kind for a matrix of size N reads 

r 2 < min{NM - 1, 1 - Nm} . 



(2.2) 



vanish. 



(2.3) 



For any stochastic matrix the characteristic polynomial is real, so we may expect a clustering of the eigenvalues of 
random stochastic matrices at the real line. This issue is related to the result of Kac, who showed that the number of 
real roots of a polynomial of order N with random real coefficients scales asymptotically like In N . The spectrum 
of a stochastic matrix is symmetric with respect to the real axis. Thus for N — 2 all eigenvalues are real and the 
support of the spectrum of the set of stochastic matrices reduces to the interval [—1,1]. 




-1 -0.5 0.5 1 




-1 -0.5 0.5 1 





FIG. 1. Eigenvalues of 3000 random stochastic matrices of size a) N=3, b) N=4. Panels c) and d) show the spectra of 
bistochastic matrices of size 3 and 4. Thick solid lines represent the bounds (2.4) of Karpelevich. 



Let Zk denote a regular polygon (with its interior) centred at with one of its fc corners at 1. The corners 
are the roots of unity of order fc. Let En represents the convex hull of En = Ufc=2-^ fc > or ' m other words, the 
polygon constructed of all fc-th roots of unity, with k = 2, . . . , N. It is not difficult to show that the support of the 
spectrum of a stochastic matrix of order N is contained in E n - see e.g. a concise proof of Schaefer |24|], p. 15 (originally 
formulated for bistochastic matrices). The fc-th roots of the unity - the corners of the regular fc-polygon, represent 
eigenvalues of non-trivial permutation matrices of size fc. For N = 3 this polygon becomes a deltoid (dotted lines in 
Fig. lb), while for N — 4 a non-regular hexagon. However, this set is larger than required: as shown by Dimitriew 
and Dynkin pq] f° r small N and later generalised and improved by Karpelevich [p7l for an arbitrary matrix size, the 
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support S^r of the spectrum of stochastic matrices forms, in general, a set which is not convex. For a recent simplified 
proof of these statements consult papers by Djokovic [gsj and Ito (29). For instance, the support Ef consists of a 
horizontal interval and the equilateral triangle (see Fig. la), while for Ef four sides of the hexagon should be replaced 
by the arcs, interpolating between the roots of the unity, and given by the solutions A = |A (t)| e^W with t £ [0, 1] of 

A 4 + (<-l)A-< = for which \(f>(t)\ £ [7r/2,27r/3], 

and 

\ 4 -2t\ 2 + (2t-t 2 -l)\{t-l)+t 2 = for which \<j)(t)\ £ [27r/3,7r]. (2.4) 



C. Spectra of bistochastic matrices 

The spectral gap of a stochastic matrix is defined as 1 — r%, where r 2 denotes the modulus of the subleading 
eigenvalue (note that in [|l2| the gap is defined as — ln^). This quantity is relevant for several applications since 
it determines the speed of the relaxation to the equilibrium of the dynamical system for which B is the transition 
matrix. Analysing the spectrum of a bistochastic matrix it is also interesting to study the distance of the closest 
eigenvalue to unity. Since for any bistochastic matrix 1 is its simple eigenvalue if and only if the matrix is irreducible, 
some information on the spectrum may be obtained by introducing a measure of irreducibility. Such a strategy was 
pursued by Fiedler fjCfl , who defined for any bistochastic matrix B the following quantities 

fi(B) := min A J2 B ^ and V ( B ) ;= mbxA J2 J2 niN-nY ^ 

where A is a proper subset of indices I = {1,2,... , N} containing n elements, 1 < n < N, and A denotes its 
complement such that A U A = I. Observe that fj, measures the minimal total weight of the 'non-diagonal' block of 
the matrix, which equals to zero for reducible matrices, while v is averaged over the number of n(N — n) elements, 
which form such a block. Fiedler (3^] used these quantities to establish the bounds for the subleading eigenvalue A2 

|1- A 2 | > 2fi(B){l -cos-^p) and |1 - A 2 | > 2v(B). (2.6) 

Since £1^ C 0^,, the supports of the spectra fulfil E^ C Ef,. Moreover, our numerical analysis suggests that for 
TV > 4 this inclusion is proper. In fact, they do not contradict an appealing conjecture [^J, that the support Ej^ of 
the spectra of bistochastic matrices is equal to the set theoretical sum Em of regular k— polygons Z^ (k = 2, . . . ,N), 
whose points are the consecutive roots of the unity. Numerical results obtained for random matrices chosen according 
to the uniform measure in the (N — l) 2 — dimensional space of bistochastic matrices are shown in Fig. lc and Id. 

It is not difficult to show that these polygons are indeed contained in the set Ej^, i.e., for each A £ [j^ =2 Zk there 
exists an N x N bistochastic matrix B such that A belongs to its spectrum, Sp(B). First note that the support _ 1 
is included in E^. Hence, it is enough to show that Zn C E^. Let us start from the following simple observation: if 
A,BG commute, then the corresponding eigenspaces are equal, and if A € Sp(A), [i e Sp(£?) are the corresponding 
eigenvalues, then x\ + (1 — x)n £ Sp(a;^4 + (1 — x)B) for each x £ [0, 1]. 

In the sequel -P(ij...ii ) will denote the matrix, which corresponds to the permutation consisting of m cycles: 

{i\ ... , (i™ —>■•••—> ), where k\ + ■ ■ ■ + k m = N. Let wn be a vector whose coordinates are the con- 

secutive TVth-roots of the unity, i.e., wn = (1, exp(27ri/iV) , . . . , exp (2ni (N — 1) /N). Let k = 0, . . . , N— 1. Then ma- 
trices P([~n) and -^(i~7V) +1 ^ comm utes, and P^~^wn — exp(2-Kki/N)wN, -P(^~/v) + ^^at = exp(27r {k + 1) i/N)wn. 
Hence the edges joining the complex numbers: exp(2irki/N) and exp(27r (k + 1) i/N) for k = 0, . . . , N — 1 are con- 
tained in E^, and so the whole boundary of the polygon Z^ . Furthermore, taking linear combinations of a bistochastic 
matrix with an eigenvalue at the edge of the polygon Zm and the identity matrix Ijv we may generate lines in E^ 
joining this boundary with point 1. It follows therefore, that the entire inner part of each polygon constitutes a part 
of the support of the spectrum of the set of bistochastic matrices. 

Let us now analyse in details the case N = 3. ^(123) denotes the 3x3 matrix representing the permutation 
1 — ► 2 — ► 3 — ► 1. Its third power is equal to identity, -P( 3 12 3) = -P(i)(2)(3) = I3 , while f( 2 12 3) = Pn.23) = -^(132)- Two 
edges of the equilateral triangle joined in unity are generated by the spectra of iP( 123 ) + (1 — x)I 3 - see Fig. 2b. 
The third vertical edge is obtained from spectra of another interpolating family iP( 123 ) + (1 — x)P^ 132 ) - see Fig. 2a. 
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The boundary of £3 is obtained from spectra of the members of the convex polyhedron of the bistochastic matrices 
of size 3. On the other hand, the permutations P(i23) and P(i2)(3) do not commute, and the spectra of their linear 
combination form a curve inside T,f - see Fig. 2c. To show that there are no points in the set T, N outside E 3 note 
that Z 2 U Z 3 = E 3 C Sf C Sf = E 3 . 

Let us move to the case N = 4. Two pairs of sides of the square forming Z4 constructed of commuting permutation 
matrices are shown in Fig. 2d and Fig. 2e, while Fig. 2f presents the spectra of a combination of noncommuting 
permutation matrices, which interpolate between 3 and 4-permutations. This illustrates the fact that E4 is contained 
in the set S^ 3 , but we have not succeeded in proving that both sets are equal. 

Analysis of the support of the spectra of stochastic and bistochastic matrices can be thus summarised by 

N 

U %k = En C Y, n C Y, n — En C En, (2.7) 

fc=2 

where En is a concave hull of the set-theoretical sum En of the regular polygons Z^ supplemented by the area 
bounded by the Karpelevich's interpolation curves |p7j-|29||, whereas En is the closed convex hull of En- 





FIG. 2. Eigenvalues of families of bistochastic matrices, which produce the boundary of the set S N . They are constructed 
of linear combinations, xP a + (1 — x)Pi,, of permutation matrices of size N = 3: a) (P(i23) and P(x32)) 5 b) (P(i23) and 73), 

c) (P(i2)(3) and h), and of size N = 4: d) (P(i 43 2) and P(i 3 )( 2 4)), e) (P(i432) and I 4 ), f) (P(i 3 42) and P( 13 2)(4))- Panels a), b), 

d) , and f) show spectra of combinations of commuting, and c) and f) of noncommuting matrices. 



III. MAJORIZATION AND ENTROPY OF BISTOCHASTIC MATRICES 



A. Majorization 

Consider two vectors x and y, consisting of N non-negative components each. Let us assume they are normalised in 
a sense that J2iLi x i = J2iLi Vi = 1- The theory of majorization introduces a partial order in the set of such vectors 
H . We say that x is majorized by y, written x -< y, if 

k k 
i=l i=l 
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for every k = 1, 2, . . . , N, where we ordered the components of each vector in the decreasing order, x± > . . . > xn 
and 2/1 > ... > y^. Vaguely speaking, the vector x is more 'mixed' than the vector, y. 

A following theorem by Hardy, Littlewood, and Polya applies the bistochastic matrices to the theory of majorization 
1- 

Theorem 1. (HLP) For any vectors x and y, with sum of their components normalised to unity, x -< y if and 
only if 

x = By for some bistochastic matrix B. (3-2) 

It was later shown by Horn [[H] that in the above theorem the word 'bistochastic' can be replaced by 'orthostochastic'. 
In general, the orthostochastic matrix B satisfying the relation x = By is not unique. 
The functions /, which preserve the majorization order: 

x <y implies f(x) < f(y). (3.3) 

are called Schur convex H. Examples of Schur convex functions include h q (x) = 5Zj=i x l f° r an y <7 > 1 , and 
= Ej=l x i ma; i- 

The degree of mixing of the vector x can be characterised by its Shannon entropy H or the generalised Renyi 
entropies H q (q > 1): 

N 

H(x) = -^Ja;* hixi, 

i=l 
1 N 

H q {x) = Y —\n{^xi). (3.4) 

1 i=l 

In the limiting case one obtains lim g _,i H q (x) = H(x). If x represents the non- negative eigenvalues of a density 
matrix, H is called the von Neumann entropy. Let x = By. Due to the Schur-convexity we have H q (x) > H q (y), since 
we changed the sign and reversed the direction of inequality. An interesting application of the theory of majorization 
in the space of density matrices representing mixed quantum states is recently provided by Nielsen jl4|]. Consider 
a mixed state p with the spectrum consisted of non-negative eigenvalues Ai > ... > Ajy, which sum to unity. This 
state may be written as a mixture of N pure states, p = Y^=\ Vi^i) if, and only if, p -< A = (Ai,... , Aat), 
where p = {pi, ... ,pj~) (if both vectors have different length, the shorter is extended by a sufficient number of extra 
components equal to zero). This statement is not true if the pure states states \ipi) are required to be distinct p2[ . 

Consider now the non-unitary dynamics of the density operators given by p — > p' = yV —1 qjUjpUj . This process, 
called random external fields p3|] , is described by L unitary operations Uj {j = 1, . . . , L) and non-negative probabilities 
satisfying Y]f—i ffj = 1- Denoting respective spectra, ordered decreasingly, by A and A' one can find unistochastic 

matrix B such that A' = BX, so due to the HLP theorem we have A' -< A (3^jlJ]. Therefore, after each iteration the 
mixed state becomes more mixed and its von Neumann entropy is non-decreasing. 

B. Preorder in the space of bistochastic matrices 

A bistochastic matrix B acting on a probability vector x makes it more mixed and increases its entropy. To settle 
which bistochastic matrices have stronger mixing properties one may introduce a relation (preordering) in the space 
of bistochastic matrices [3f| writing 

B\ < £>2 iff B\ = BBi for some bistochastic matrix B . (3-5) 

We have already distinguished some bistochastic matrices: permutation matrices P with only N non-zero elements, 
and the uniform van der Waerden matrix B*, with all its N 2 elements equal to 1/N. For an arbitrary bistochastic 
matrix B and for all permutation matrices P we have £?* = B*B and B = (i?P _1 ) P, and hence B* < B < P . The 
relation B\ -< B2 implies B\y -< B2y for every probability vector y, but the converse is not true in general (see |3(| 
and 0). 
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C. Entropy of bistochastic matrices 



To measure mixing properties of a bistochastic matrix B of size N one may consider its entropy. We define Shannon 
entropy of B as the mean entropy of its columns (rows), which is equivalent to 



N N 

//::/>>:: y X 1:1 ^ • ( 3 - 6 ) 

1=1 j=l 

As usual in the definitions of entropy we set OlnO = 0, if necessary. The entropy changes from zero for permutation 
matrices to In TV for the uniform matrix B* . 

Due to the majorization of each column vector, the relation C < B implies H(B) < H(C), but the converse is not 
true. To characterise quantitatively the effect of entropy increase under the action of a bistochastic matrix B, let us 
define the weighted entropy of matrix B with respect to a probability vector y = (yi, . . . , i/n)- 



N N N 

fc=l k=l j=l 

where B^ is a probability vector defined by Bk := (-Bife, ■ ■ ■ , B^k) for fe = 1, . . . , N . In this notation H(B) = H Ct (B), 
where e* = (1/N, . . . , 1/N). The weighted entropy allows one to write down the following bounds for H(By) 



max {H(y),Htf(B)} < H(Bij) < H(y) + H ff (B) . (3.8) 

the proof of which is provided elsewhere ]3§| ]. These bounds have certain implications in quantum mechanics. For 
example, if a non-unitary evolution of the density operator under the action of random external field is considered 
fl4| , they tell us how much the von Neumann entropy of the mixed state may grow during each iteration. 

Using the above proposition we shall show that the entropy of bistochastic matrices is subadditive. Namely, the 
following theorem holds: 

Theorem 2. Let A and B be two bistochastic matrices. Then 



max {H \A), H(B)} < H(AB) < H{A) + H{B) , (3.9) 

and, analogously, 

max{H(A),H(B)} < H(BA) < H(A) + H(B) . (3.10) 

Proof. We put C := AB a nd consider stochastic vectors y n := (ci n , . . . , cjv n ) for n — 1, . . . ,N (the columns of 
the matrix C). Applying (|3.8|) we get 



max {Hff n (B),H {y n )) < H (By n ) < H ffn (B) + H (y n ) . 

Hence 

max (j2k=l c knJ2f=lV {bjk) , Efcl 7 ? (Cfcn)) < J2k=lV [EjLlhjCjn 

and 

J2k=lV (X^LAjCj") < J2k=l C knJ2f=lV(bjk) +Ef=i 7 ?( c fcn) , 

wher e n (x) = — xlnx for x > 0. Mult iplying the above equalities by 1/N and summing over n = 1, . . . , N we get 
(3.9). Setting C — BA we obtain ( 3.10| ) in an analogous way. □ 
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IV. ENSEMBLES OF RANDOM UNISTOCHASTIC MATRICES 



A. Unistochastic matrices 



To demonstrate that a given bistochastic matrix B is unistochastic one needs to find unitary matrix U such that 
, | 2 . In other words one needs to find a solution of the coupled set of nonlinear equations enforced by the 
unitarity UW = WU = I, 



Bij - 



N 

E 



y/B^Bji exp(i0 jfe - 4>j{] = S k i for all 1 < k < I < N 



(4.1) 



for the unknown phases of each complex element Uj k — 
since B is bistochastic. 



The diagonal constraints for k = I are fulfilled, 



n%. This is not the 



For N = 2 all bistochastic matrices are orthostochastic (see eg. Eq.(4.6)), and so = fi^ 
case for higher dimensions, for which there exist bistochastic matrices which are not unistochastic, and so f2 
for N > 3. Thus is not a convex set for N > 3, since it contains all the permutation matrices and is smaller 
than their convex hull, which, according to the Birkhoff theorem, is equal to f2|J . 



N 



matrices which are not unistochastic were already provided (for N 



Simple examples of bistochastic 
3) by Schur and Hoffman 



B x 



1 1 
1 1 
1 1 



B, 



3 3 
3 1 2 
3 2 1 



(4.2) 



b) 




FIG. 3. Chain rule for unistochasticity for N = 3: (see (4. Z) and (i.4)), a) the longest link L\ > Lz + L3 so the matrix B 
is not unistochastic, b) condition for orthostochasticity L\ — L2 + I/3, c) weaker condition for unistochasticity L\ < L2 + L3. 



To see that there exists no corresponding unitary matrix U we shall analyse constraints implied by the unitarity. 
Define vectors containing square roots of the column (row) entries, e.g., ilk := {\/Bki, VBk2, V-Sfejv}- The scalar 
products of any pair of any two different vectors Vk ■ vi consists of N terms, L„ = y/Bk n Bi n , n = 1, ...,N. In the 
case of B\ from (4.2) the scalar product related to the two first columns (vi,V2) consists of three terms L\ = 1, 
L2 = L3 = 0, which do not satisfy the triangle inequality. Thus it is not possible to find the corresponding phases 
(f>jk satisfyig (4.1). This observation allows us to obtain a set of necessary conditions, a bistochastic B must satisfy 
to be unistochastic Jill : 



max ^ \J B mk B m i < - ^ y/B jk Bji , for all 1 < k < I < 



N 



(4.3) 



and 



1 N 

max ^ y/B km Bim < - y/B kj Bij , for all 1 < k < I < 



N . 



(4.4) 



We shall refer to the above inequalities as to a 'chain rule': the longest link L\ of a closed cha in cann ot b e longer 
than the sum of all other links L2 + ... + Ln - see Fig. |^. The set of N(N — l)/2 conditions (4.3) (resp. (4.4)) treats 
all possible pairs of the columns (resp. rows) of B. Only for TV = 3 the both sets of conditions are equivalent (since 
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the equations (4.1) for the phases may be separated |p9| ), but for N > 4 there exist matrices which satisfy only one 
class of the constraints. For example, a bistochastic matrix 





" 1 


17 


25 


57" 


1 


38 


38 


24 





100 


42 


5 


29 


24 




. 19 


40 


22 


19 



(4.5) 



found by Pakohski g^l, satisfies the row conditions (4.4), but violates the column conditions (4.3), so cannot be 
unistochastic (the third term of the product of the roots of the first and the fourth column is larger than the sum 
of the remaining three terms). It is easy to see that for an arbitrary TV these necessary conditions are satisfied by 
any bistochastic B sufficiently close to the van der Waerden matrix for which all links of the chain are equal, 
Li — 1/N. It is then tempting to expect that there exists an open vicinity of B* included in f2^, i.e., that -B* lies in 
the interior of £1^ and, consequently, that il K has positive volume. Certain conditions sufficient for unistochasticity 
were already found by Au-Yeung and Cheng J40|, but they do not answer the question concerning the volume of f2^. 
On the other hand, it is well known that the set of unistochastic matrices is connected and compact pjj , and is not 
dense in the set of bistochastic matrices [H . 

To analyse properties of bistochastic matrices it is convenient to introduce so called T-transforms, which in a sense 
reduce the problem to two dimensions. The T-transform acts as the identity in all but two dimensions, in which it 
has a common form of an orthostochastic matrix 



cos 2 <p sin 2 tp 

• 2 9 

sin ip cos ip 



such that O2 (<p) — 



cos if sin ip 
— sin <p cos ip 



(4.6) 



is orthogonal. Any matrix B obtained as a sequence of at most (JV — 1) T— transforms, B = Tjv_i • ■ - T2T±, where 
each Tk acts in the two-dimensional subspace spanned by the base vectors k and k + 1, is orthostochastic. To show 
this it is enough to observe that each element of B is a product of non-trivial elements of the transformations Tfc. 
Hence taking its square root and adjusting the signs one may find a corresponding orthogonal matrix defined by the 
products of the elements of Ok Although any product of an arbitrary number of T-transforms satisfies the 

chain-links conditions |32| , for TV > 4 there exist products of a finite number of T-transforms (also called pinching 
matrices) which are not unistochastic p3| ] . In the same paper it is also shown that there exist unistochastic matrices 
which cannot be written as a product of T-transforms. 

Consider a unistochastic matrix B and the set Mb C U(N) of all unitary matrices corresponding to B in the sense 
that Bij = \Uij\ 2 for i,j = 1, ... , N. Such sets endowed with appropriate probability measures play a role in the 
theory of quantum graphs [^0|-^2| and were called unitary stochastic ensembles by Tanner |0] . It is easy to see that 
these sets are invariant under the operations U — * V1UV2, where V\ and V% are arbitrary diagonal unitary matrices. 
The dimensionality arguments suggest that, having U € Ub fixed, each element of Ub can be obtained in this way 
p4j . However, in general this conjecture is false and for certain bistochastic matrices B the set Ub is larger. This is 
shown in Appendix |a|, in which a counterexample for N = 4 is provided. 



B. Definition of ensembles 



To analyse random unistochastic matrices one needs to specify a probability measure in this set. As it was already 
discussed in the introduction, unistochastic (USE) and orthostochastic (OSE) ensembles can be defined with help of 
the Haar measure on the group of unitary matrices U(N) and orthogonal matrices O(N), respectively [^9| . In other 
words the bistochastic matrices 



:= l^f, and B% := (O^) 2 , (4.7) 

pertain to USE and OSE respectively, if the matrices U and O are generated with respect to the Haar measures on 
the unitary (orthogonal) group. 

Dynamical systems with time reversal symmetry are described by unitary symmetric matrices [ f45| . The ensemble 
of these matrices, defined by W — UU T , is invariant with respect to orthogonal rotations, and is called circular 
orthogonal ensemble (COE). In an analogous way we may thus define the following three ensembles of symmetric 
bistochastic matrices (SBM) 
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a) 5! 

b) S 2 

c) S 3 



= BB T ; so Sij :=Efc=il^fc| 2 |;7 jfe | 2 , 

= \(B + B T ); so Sa = ±(\Uij\ 2 + |l^| 2 ), 

= \w tJ \ 2 = \(uu T ) v \ 2 ; so S^ :=\EtiU ik U jk \ 2 , 



(4.8) 



where bistochastic matrices B are generated according to USE (or, equivalently, unitary matrices U are generated 
according to CUE). 



C. Spectra of random unistochastic matrices 

Since the sets of the bi-, uni-, and (ortho-) stochastic matrices are related by the inclusion relations: 
C C n% C tiff, analo gous relations must hold for the supports of their spectra, C C C E^-. 
For N — 2 the spectrum of a bistochastic matrix must be real. The subleading eigenvalue A2 = 2Bu — 1, which 
allows one to obtain the distributions along the real axis. For USE the respective density is constant, P r (X) = 1/4 
for A G (—1, 1), while for OSE it is given by the formula, P r (\) = 1/(2tt\/I — A 2 ). These densities are normalized to 
1/2, since for any random matrix its leading eigenvalue is equal to unity. 

For larger matrices we generated random unitary (orthogonal) matrices with respect to the Haar measure by means 
of the Hurwitz parameterisation |46| as described in W7p8]. Squaring each element of the matrices generated in 
this way we get random matrices typical for USE (OSE). In general, the density of the spectrum of random uni-, 
(ortho-)stochastic matrices may be split into three components: 

a) two-dimensional density with the domain forming a subset (Ejv) of the unit circle; 

b) one-dimensional density at the real line described by the function P r (x) for x G [—1,1], 
and 

c) the Dirac delta jjS(z — 1) describing the leading eigenvalue, which exists for every unistochastic matrix. 





-1 -0.5 0.5 1 



-1 -0.5 0.5 1 





-1 -0.5 0.5 1 



-1 -0.5 0.5 1 



FIG. 4. Spectra of random unistochastic matrices of size a) N = 3 (1200 matrices) and b) N = 4 (800 matrices); spectra 
of random orthostochastic matrices of size c) TV = 3 (1200 matrices) and d) N = 4 (800 matrices). Thin lines denote 3- and 
4-hypocycloids, while the thick lines represent the 3-4 interpolation arc (part of it is shown in Fig. 2f). 
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Basing on the numerical results presented in Fig. 4a, we conjecture that £3 = £3 and consists of a real interval 
(already present for N = 2) and the inner part of the 3-hypocycloid. This curve is drawn by a point at a circle of 
radius 1/3, which rolls (without sliding) inside the circle of radius 1. The parametric formula reads 



= -3- (2 cos < 



cos 20), 
y = i (2 sin — sin 20), 



(4.9) 



where <j) € [0, 2n). 

To find the unistochastic matrices with spectra at the cycloid consider a two-parameter family of combinations of 
permutation matrices a 2 I + b 2 P + c 2 P 2 with a 2 + b 2 + c 2 = 1. Here P represents the nontrivial 3 elements permutation 
matrix, P — P(i23), so P 3 = I. One-parameter family of these bistochastic matrices, which satisfy the condition for 
unistochasticity, produce spectra located along the entire hypocycloid. Consider the matrices 



3 (<p) ■■= 



a b c 
cab 
b c a 



and B^{ip) := 



b 2 c 2 

a 2 b 2 

9 9 

c a z 



(4.10) 



where their elements depend on a single parameter ip £ [0, 2tt) and 



'a = a((f) = — 1(1 + 2 cosy?) 

b= %) = i(cosy5-l) + -^=siny? , (4.11) 
c(ip) = i(cosy3 - 1) - 4= siny? 



It is easy to see that a 2 + b 2 + c 2 = 1 and ab + be + ca — 0, so O3 is orthog onal and B3 is orthostochastic. Simple 
algebra shows that the spectrum of B^(ip) forms the 3-hypocycloid given by (4.9) - see Appendix |b| 



An alternative approach, based on exponentiation of permutation matrices, leads to unistochastic matrices with 
spectrum on a hypocycloid. Let Pjq :— P(i2—n) be the nontrivial permutation matrix of size N. Since Pjy is unitary, 
so is its power (P/v) Q . We define it for an arbitrary real exponent by (-P/v) a := W D a U, where U is an unitary matrix 
diagonalizing P/v and D represents the diagonal matrix of its eigenvalues - it is assumed that the eigenphases of such 
a matrix belong to [0, 27r). Since Pjy = Ij\r, one may expect that defining the corresponding bistochastic matrices 
and varying the exponent a from zero to unity one obtains an arc of a hypocycloid. This fact is true as shown in 
Appendix |C|, in which we prove a general result, valid for an arbitrary matrix size. 

Proposition 3. The spectra of the family of unistochastic matrices of size N 



C 1 H(^)«r with «e[0,A/2], (4.12) 
generate the N -hypocycloid (and inner diagonal hypocycloids of the radius ratio k/N with k = 2, . . . , N — 1). 




FIG. 5. Spectra of a family of bistochastic matrices B^ N ' a ' with a € [0, N/2] obtained by exponentiation of the permutation 
matrix Pjv plotted for a) N = 5 (hypocycloids 5 and 5/2), b) N = 6 (hypocycloids 6, 3 = 6/2 and 2 = 6/3), and c) N — 7 
(hypocycloids 7, 7/2 and 7/3). 
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Let Hn denote the set bounded by N— hypocycloid. Fig. 4c suggests that P3 is contained in E^. This conjecture 
may be supported by considering the spectra interpolating between the origin (0, 0) and a selected point on the hypocy- 
cloid. To find such a family we shall use the Fourier matrix F^ N ^ of size N with elements := exp(— 2klwi/ N) / y/N . 
Since amplitudes of all elements of this matrix are equal, the corresponding bistochastic matrix equals to the van 
der Waerden matrix for which all sublcading eigenvalues vanish. The matrix Pjy generates a unistochastic matrix 

B( N > a ) with an eigenvalue at the hypocycloid, so the family of unistochastic matrices related to (P^)^ '(F^) 1 13 
provides an interpolation between the origin and a selected point on the hypocycloid - see Fig. 6. In other words, we 
conjecture that the spectra of a two parameter family of unistochastic matrices obtained from (Pjv)" (F^ N ^) b explore 
the entire Hjq- 

Numerical results obtained for N — 3 random uni-, (ortho-)stochastic matrices allow us to claim that there are no 
complex eigenvalues outside the 3- hypocycloid, so that E^ = = P2 U H 3 . Interestingly, H 3 - the 3- hypocycloid 
and its interior, determines the set of all unistochastic matrices which belong to the cross section of Cl§ defined by 
the plane spanned by P3 , Pf and P| = I |32j , see also Appendix |b[ 




FIG. 6. Spectra of a family of bistochastic matrices obtained be squaring elements of unitary matrices f-Pj?) (F^-^) 1 13 
with fl £ [0, 1] and a) N = 3, the curves are labelled by a = 0, 1/8, 2/8, .. . , 3/2, b) N = 4, a = 0, 1/8, 2/8, . . . ,2. For reference 
we plotted the hypocycloids with a thin line. 



Numerical results for spectra of random uni-, (ortho-)stochastic matrices for N = 4 are shown in Fig. 4b and 4d. 
The set Ejjf includes the entire set Eg 7 , but also the 4-hypocycloid, sometimes called asteroid, and formed by the 
solutions of the equation 

;E 2/3 +y 2/3 = 1 (413) 

More generally, the spectra of orthostochastic matrices of size N contain the iV-hypocycloid Hm 



X = i[(JV-l)cOB0 + COB(JV-l)0], 

y = i[(iV-l)sin0-sin(7V-l)0], 

with corner at 1, as it is proved in Appendix [B|. The set E^ contains of the sets bounded by hypocycloids of a smaller 

dimension, which we denote by Gn :— Ufc!=2 ^k- 

However, as seen in Fig. 4b there exist some eigenvalues of unistochastic or even orthostochastic matrices outside 
this set. In fact one needs to find interpolations between roots of unity of different orders (the corners of k and 
n-hypocycloid) based on the families of orthostochastic matrices. Such a family interpolating between the corners of 
H3 and H4 is plotted in Fig. 2f, since these bistochastic matrices are orthostochastic. A general scheme of finding the 
required interpolations is based on the fact that, if 0(cp) is a family orthogonal matrices and Bij(ip) = O^(ip) is the 
corresponding family of orthostochastic matrices, than another such family is produced by the multiplication by some 
permutation matrix P: 0'(tp) = 0(tp)P. F or ex ample, if the matrix Oi(ip) of size 4 contains elements On = O22 = 1 



and the block diagonal matrix 02(f) (see (4.6)), then the squared elements of the orthogonal matrices 04(y)P(i234) 
provide orthostochastic matrices, the spectra of which form the 3-4 interpolating curve contained in E4 and plotted in 
Fig. 2f and 4d. To obtain orthostochastic matrices, the spectra of which provide the curve which joins both complex 
corners of H 3 and does not belong to it (see Fig 4d), one should take the orthogonal matrix 31 of size 4 with On = 1 



which contains the block diagonal matrix Oa((p) (4.1C), and create bistochastic matrices out of 03.i(v)P(i2)(34)- The 
exact form of these families of orthostochastic matrices is provided in Appendix |)| 
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Let Gm denote the set Gn extended by adding the regions bounded by the interpolations between all neighbouring 
corners of Hk and H n with k < n < N constructed in an analogous way as for N = 4. The sequence of hypocycloids 
with the neighbouring corners is given by the order of fractions present in the Farey sequence and Farey tree jl9| . 
Our investigation of t he s upport of the spectra of unistochastic and orthostochastic matrices may be concluded by a 
relation analogous to (EJ) 



N N 

(J H k := Gn C S° ~ S^v — &n C Bjv =: |J ^fe, (4.15) 

k=2 k=2 

where the sign ~ represents the conjectured equality. Numerical results suggest that the support Sat is the same for 
both ensembles: USE and OSE. However, the density of eigenvalues P(z) is different for uni- and (ortho-)stochastic 
ensembles. The repulsion of the complex eigenvalues from the real line is more pronounced for the orthostochastic 
ensemble, as demonstrated in Fig 4b and 4d. 

The larger N, the better the domains T, N and T, N fills the unit disk. On the other hand, the eigenvalues of a large 
modulus are unlikely; the density P(z) is concentrated in a close vicinity of the origin, z = 0. Moreover, in this case 
the weight of the singular part of the density at the real line, decreases with the matrix size N. To characterise the 
spectrum qualitatively we analysed the densities of the distributions Pk if) of the moduli of the largest eigenvalues Xk ■ 
The results obtained for random unistochastic matrices are shown in Fig. 7a. Fig. 7b shows analogous data for the 
singular values er,; of B - per definition square roots of the real eigenvalues of the symmetric matrices BB T [^o| . Since 
the singular values bound moduli of eigenvalues from above [ pCf l, the distribution P{cr) is localised at larger values 
than P(r). Thus the expectation values satisfy (rfe) < (07.) as shown in Fig. 7c and 7d. The modulus of the second 
eigenvalue decreases with matrix size as N^ 1 / 2 . These results are consistent with the recent work of Berkolaiko [ fl9| |, 
who suggested describing the distribution Pzir) by the generalised extreme value distribution pi)] . 




FIG. 7. Unistochastic ensemble. Distribution of the modulus of the second (A), third (□), fourth (O) and fifth (x) 
largest eigenvalues a), and singular value b) (lines are drawn to guide the eye). Expectation value (r 2 ) for k = 2, 3, 4, 5 of the 
subleading c) eigenvalues, and d) singular values as a function of the matrix size N. Exponential fit, represented by solid lines, 
give (r 2 ) w exp(-0.503) and (<r 2 ) ~ exp(-0.446) 
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D. Average entropy 



We can compute the mean entropy averaged over the ensembles of unistochastic or orthostochastic matrices. Clearly 
we have 

(H)use =N-{- |t/ 4 /m|^| 2 ) (4.16) 

and 



(H) OSE = N-(-\O ij \ 2 ln\O ij \ 2 ) , (4.17) 

where Uij and Oij are entries of unitary or orthogonal matrices, respectively. The exact formulae for the above 
averages were obtained by Jones in |5^,^3| . Using these results we get 

N 1 

(tfW = *(^+i)-*(2) = Ei ( 4 - 18 ) 



k 

k=1 



and 



( flW ^WH.)-W)^^-' N N Zf h + 1 ■ <«») 

where "J denotes the digamma function. For large N the digamma function behaves logarithmically, ^f(N) ~ In TV, 
and both averages display similar asymptotic behaviour (H) USE « In N — 1 + 7 and (H) OSE w In N — 2 + 7 + In 2 , 
where 7 w 0.577 stands for the Euler constant. Note that both averages are close to the maximal value In A, attained 
for orthostochastic matrices corresponding to and their difference (H) USE — (H) OSE converges to 1 — In 2 0.307. 



E. Average traces 

Traces of consecutive powers of bistochastic matrices, tiB n , are quantities important in applications related to 
quantum chaos ||. In the following we evaluate the traces tjv, n := (tr£?")usE, averaged over the unistochastic 
ensemble. 

The spectrum of a bistochastic matrix size 2 is real and may be written as Sp(-B) = {l,y} with y G [—1, 1]. For 
the unistochastic ensemble y — cos 26, where the angle 6 is distributed uniformly, P{9) = 2/tt for 9 £ [0, 7r/2]. The 
average traces may be readily expressed by the moments (y k ) for k £ N, namely, 

t 2 ,2k+i = 1 + (y 2k+1 ) = 1 , 
h, 2k = 1 + (y 2k ) = §±f . (4.20) 

Since the average size of all the diagonal elements must be equal, (\Un\ 2 ) = l/N, so the average trace £jvi = 
(J2iLi l^iil 2 ) = 1 an d does not depend on the matrix size. Using the results of Mello 54|, who computed several 



averages over the Haar measure on the unitary group U (N) , we derive in Appendix [|] the following formulae 

1 , f 1 for N = 2 , , A . 

tN > 2 = l + —y and ^ 3 = \i + i^ot fOT ^^ 3 - ( } 

Analogous expressions for larger N may be explicitly written down as functions of the Mello averages, but it is not 
simple to put them into a transparent form. Numerical results support a conjecture that for arbitrary N the average 
traces tend fast to unity and the difference ijv,n — 1 behaves as A 1 ^™. This fact is related to the properties of 
the spectra discussed above: for large N the spectrum is concentrated close to the center of the unit circle, so the 
contribution of all the subleading eigenvalues to the traces becomes negligible. Related results on average traces of 
the symmetric unistochastic matrices BB T are provided in |19| . 
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V. CLOSING REMARKS 



In this paper wc define the entropy for bistochastic matrices and proved its subadditivity. We also analysed special 
classes of bistochastic matrices: the unistochastic and orthostochastic matrices, and introduced probability measures 
on these sets. We found a characterization of complex spectra of unistochastic matrices in the unit circle and discussed 
the size of their subleading eigenvalue, which determines the speed of the decay of correlation in the dynamical system 
described by the matrices under consideration. Unistochastic matrices find diverse applications in different branches 
of physics and it is legitimate to ask, how a physical system behaves, if it is described by a random unistochastic 
matrix Since we succeeded in computation of mean values of some quantities (traces, entropy) averaged over the 
ensemble of unistochastic matrices, our results provide some information concerning this issue. 

On the other hand, several problems concerning bi-, uni-, and (ortho-)stochastic matrices remain open and we 
conclude this paper listing some of them: 

a) Prove the conjecture that the union of the spectra of bistochastic matrices E^ is equal to the sum of regular 
polygons E N , i.e., E^ = E N . 

b) Prove that the support Eg' of spectra of unistochastic matrices of size N = 3 contains exactly the interval [—1,1] 
and the set bounded by the 3-hypocycloid H 3 , i.e., E^ = G 3 . 

c) Find an analytical form of the boundaries of the set Gn obtained of interpolations between the neighbouring 
corners of fc-hypocycloids with k < N . 

d) Show that the support E^ of the spectra of unistochastic matrices contains Gn and check whether both sets 
are equal, i.e., E^ = Gn- 

e) Show that the supports of spectra of uni- and (ortho-)stochastic matrices are the same, i.e., E^ = E^y. 

f) Calculate probability distribution Pn{z) of complex eigenvalues ensembles for uni-(ortho-)stochastic matrices. 

g) Compute the expectation value of the subleading eigenvalue (r 2 ), averaged over uni-(ortho-)stochastic ensemble. 

h) Calculate averages over the unistochastic ensemble of other quantities characterizing ergodicity of stochastic 
matrices, including ergodicity coefficients analysed by Seneta J55| and entropy contraction coefficient introduced by 
Cohen et al. Q. 

i) Find necessary and sufficient conditions for a bistochastic matrix to be unistochastic. 

j) Provide a full characterization of the set of all unitary matrices U, which lead to the same unistochastic matrix 
B, i.e., Bij = \Uij\ 2 for i, j = 1, . . . , k. 

This paper is devoted to the memory of late Marcin Pozniak, with whom we enjoyed numerous fruitful discussions on 
the properties of bistochastic matrices several years ago. We are thankful to Prot Pakohski for a fruitful interaction 
and also acknowledge helpful remarks of I. Bengtsson, G. Berkolaiko, G. Tanner, and M. Wojtkowski. Financial 
support by Komitet Badah Naukowych under the grant 2P03B-072 19 and the Sonderforschungsbereich 'Unordung 
und grosse Fluktuationen' der Deutschen Forschungsgemeinschaft is gratefully acknowledged. 



APPENDIX A: UNISTOCHASTIC MATRICES STEMMING FROM A GIVEN BISTOCHASTIC MATRIX 



Let us consider two unitary N x N matrices U and W such that for all i, j = 1, . . . , N, 



|t%| 2 = l^f , (Al) 



i.e., that the corresponding unistochastic matrices are the same. It is obvious that this happens if W = V\UV^ 



2 



with V\ and V<i unitary diagonal. However the converse statement, i.e., that (Al) implies existence of two diagonal 
unitary matrices V\ and Vi such that U = V1WV2, is false p4j . The plausibility of such conjecture is based on the 
following dimensional argument: we have TV 2 real numbers := |fZy| fulfilling 2N — 1 relations stemming from the 
normalisation of the rows: 



5>« = 1. i = l,...,N, (A2) 

and the columns 
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N 



= 1, 3 = 1,. 



(A3) 



of th e un itary matrix U. The number of in dep endent relations is less by one than the total number of eq uations in ( A2) 
and ( |A3| ) since summing all equations in ( A2) over i gives the same as summing all equations in (A3) over j, namely 
N = Tr(U^U). On the other hand the left and right multiplications by unitary diagonal V\ and V% introduce exactly 
2 A — 1 parameters (here the number of the independent parameters is diminished by one from the number of non-zero 
elements of both D\ and D 2 since in the resulting matrix only the differences of eigenphases of D\ and D% appear, 
so we can always put one of the eigenphases of D±, say, to zero without changing the result of the transformation 
W i— > V1WV2). The simplest counterexample we know involves the following unitary matrices U and W 





" 1 


1 


1 


1 " 


1 


1 


-1 


-e la 


e la 


2 


1 


-1 


e la 


-e la 




. 1 


1 


-1 


-1 



w = ~ 



1 

-1 



1 

-1 

,03 
oi/3 



(A4) 



with \Ui, 



m. 



1/4. It is a matter of simple explicit calculations to show that there are no unitary diagonal 



Vi and V 2 fulfilling W = V1UV2, if a, /3 G [0, 2tt] and af3 ^ 0. 



APPENDIX B: ORTHOSTOCHASTIC MATRICES WITH SPECTRUM AT HYPOCYCLOIDS 

In this appendix we construct orthostochastic matrices with spectra on A^-hypocycloids. Consider the following 
N x N permutation matrix 



P : = 



1 






1 








1 




(Bl) 



where for simplicity the dimensionality index N has been omitted. We have P N = I, P K 7^ I for K < N and the 
eigenvalues of P equal exp (^fc), k = 0, 1, . . . , N — 1. 

We shall discuss separately the cases of odd and even N. 

First let N — 2K + 1 and 



JV-l JV-l 

3=0 3=0 



(B2) 



Observe that, since P m and P n do not have common non-zero entries for m 7^ n, < to, n < N — 1, the elements of 
B are squares of the corresponding elements of O. The eigenvalues of O and B read, respectively, 



N-l 



3=0 
N-l 



2m 



Ak = £ °J exp ( "/v"^ ) ' = 0, 1, . ■ • , - 1 



Afe = J2 a ) ex p ( — kj ) ' * = 0, i, . . . , - 1. 

3=0 



Inverting the discrete Fourier transforms in (B3) we obtain 



N-l 



= ^E Afeex p(-^^)' k = 0,1,.. . , A — 1 



fc=0 



(B3) 
(B4) 



(B5) 



1G 



which, upon substituting to (B4), gives the eigenvalues of B in terms of the eigenvalues of O 

N-1N-1N-1 / g ^ v ^ JV-liV-1 ^ N-l 

Xk =mY,Y, E AzA r exp(-^(Z + r-fc)j) = ^E E A * A A,fc-i = ^ E A ' A *-«' ( B6 ) 



j=0 1=0 r=0 



1=0 r=0 



1=0 



where the indices are counted modulo N. 

Our aim is now to find a family of orthogonal matrices 0(<p) such that when 4> changes (from to 2ir, say) the 
eigenvalue Aq(</>) renders the .ZV-hypocycloid Hn in the complex plane, i.e. 



First observe that the desired result is achieved if 



(B7) 



Indeed 



A k = exp[i(k - K)<f>], k = 0, 1, . . . , 2K = N - 1. 



JV-l 



JV-l 



N-l 



A ° = ^ E AA -< = A o + E AA - = A o + E AA --< 

i=0 \ J=l / \ 1=1 / 

= J_ ( e - 2 ^ + ^ e i(l-K)<j> e liN-l-K)<l> j = J_ ^ e -i(JV-l)^ + j-jy _ ^ e i(JV-2A-)^ 



(iV - l)e 40 + e"*^- 1 ^ 



(B8) 



(B9) 



(BIO) 



In order to construct an orthogonal matrix O with the spectrum (B8) it is enough to find an antisymmetric A with 
the eigenvalues 



H k = i(k-K), k = 0,l,...,2K. 



(Bll) 



If A is a polynomial in P then so is O = exp(A</>), moreover O is orthogonal and has the desired spectrum (B8). Let 
us thus write 



•1: E^" /,A 

3=1 



which is clearly antisymmetric, with the eigenvalues 



(B12) 



K 
3=1 



. 2ni ., \ / 2m 

exp I — ]k I -exp I — (N - j)k 



K 

= ^E 

3 = 1 



2tt 

a 9 - sin , 
J y 2K + 1 



kj), k — 0,1, . . . , 2K. (B13) 



Obviously /iq = and fi^K+l-k = fJ-k, k = 1,2, . . . , K. To fulfil (Bll) it is thus enough that 



K 



2 > a ,- sin , 

3 = 1 



2tt 



fcj =fc, fc = 1,2,... 



(B14) 



Using 



2tt 



A" 

E si n„>A - 1 

3 = 1 



kj I sin 



2tt 



2if+l 



2K + 1 



(B15) 
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we solve (B14) for ay, 



K 



2K 



iE fcsi 



k=l 



2ir . 
2K + 1 J 



(B16) 



For N even, TV = 2K, the construction is very similar. We introduce E := diag(l, 1, . . . ,1, —1) and define P := EP 
which has eigenvalues exp {^ff{k + ^)) , k = 0, 1, . . . , N — 1, and construct the matrix O as a polynomial in P rather 
than in P 



with the eigenvalues: 



N-l 

0:=E "•• / "- 

i=o 



(B17) 



Af-1 



A fe = J2 a i ex P (^T ( fc + 5 ) J ) > fc = 0, 1, . . . , JV - 1 



The matrix B of the squared elements of O is, as previously, the following polynomial in P, 

N-i 

Using exactly the same method as above we express the eigenvalues of B are given in terms of Aj 

1 Ar ~ 1 . . 
z=o 

Now if 



(B18) 



(B19) 



(B20) 



A fe = exp 



i \ k — K + - 



, fc = 0,l,... ,2K-1 = N-1, 



then 



-i(2AT-l)0 



JV-1 

E< 



(B21) 



+1/2)0 (N-l-K+\/2)4> 



= 1 ( e -i(iv-i)0 + {N _ 1)e i(N-2K+i)^ = I [(AT - l) e ** + e -<(^-i)*] . 
In full analogy with the case of odd N we look for an antisymmetric matrix A with the eigenvalues 



(B22) 



(Et fe = i [k - K 



(B23) 



in the form of a polynomial in P 



A := 2 1 / 2 5kP A ' + J2 & i { P + pN ~ 3 ) 



(B24) 
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Since, as it is easy to check, (P-') T = —P N J for j = 1, 2, . . . , K, the matrix A is indeed antisymmetric for arbitrary 
real ay, j = 1, 2, . . . , K. The eigenvalues of A read: 



K-l 

(ik = 2 1/2 i(-l) k a K + 5 J sin 

3=1 







[s(* H 







fc = o,i,... ,2/r-i. 



Using arguments similar to those in the odd N case, we conclude that the choice 

K-l 



ft, 



a K 



— sir 
K ^ 

4 ' 



leads to the desired result ( B23 ) and, consequently, the 2JC-hypocycloid ( |B22| ). 



(B25) 



(B26) 
(B27) 



APPENDIX C: UNISTOCHASTIC MATRICES WITH SPECTRUM AT HYPOCYCLOIDS 

The above constructed matrices with spectra on JV-hypocycloids were ortho 'stochastic. If the desired matrix should 
be merely unistochastic, but not necessarily orthostochastic, the construction is even simpler. To this end let us 
consider the matrix 



P a := U^D a U, 



(CI) 



where U is an unitary matrix diagonalizing P, where P is given by (Bl)), and D is a diagonal matrix 



D := diag (l, e 2 ^ N , e 4 "/ N , . . . , e ^-i)*i/^ 
with the eigenvalues of P on the main diagonal. Hence, consequently 



(C2) 



D a =diag(A ,Ai,... ,A N _i), 



(C3) 



where A& := exp(2/ca7ri/7V), k = 0, 1, . . . , N— 1, are the eigenvalues of P a . Obviously P a is unitary, and as a function 



of P can be written in the form P a := X^o* a i^ >3 ' > which gives for the eigenvalues 



N-i 



A k = e «-i/iv = £ a . exp / ™ kj \ 5 fc = 0, 1, . . . ,N — 1 



N 



(C4) 



As previously we obtain the coefficients aj by inverting the discrete Fourier transform 



N-l 



1 C / 2ni . 

2^A fe exp( — —kj), A: = 0,1,. 



N 



k=0 



N 



,N-1. 



(C5) 



The associated unistochastic matrix reads thus 



N-l 

3=0 



(C6) 
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and has as the eigenvalues 



N-l 



' exp 



2m, . 

hi 

N J 



k = 0,1,... ,N- 1, 



(C7) 



i.e. 



j=0 1=0 r=0 



N-l N-l 



N-l 



N-l N-l N-l / O • \ 1 I 

^ = iEEE A <^ «P + r- k)j = 1 £ £ W r ,_ fc = i £ AjA*_ 



1=0 r=0 



(C8) 



hence 



Ai 



A A 



JV-l 



E a ' a ^-a 



i=l 



1 

iV 



/ e -2(jv-i)Q7ri/JV + (TV - i) e 2 <»™/^ 



(C9) 



which renders the iV-hypocycloid for a G [0, N/2]. Similarly, the further eigenvalues Xk generate the inner hypocycloids 
(e.g. 3-hypocycloid for N = 6 - see Fig. 5), what proves Proposition 3. 

Observe that for N — 3 the permutation matrices P, P 2 and P 3 = I3 form an equilateral triangle (in sense of the 
Hilbcrt-Schmidt distance, which is induced by the Frobenius norm of a matrix, ||P|| := \/PPt). Comparing Eq. (C6) 
and (C7) with k = 1 we see that both quantities have the same structure and the same dependence on the coefficients 
Oj, which are implicit functions of the parameter a. Therefore, varying this parameter we obtain the very same curves 
in two entirely different spaces: the eigenvalue Ai = Xi(oe) provides the 3-hypocycloid in the plane of complex spectra, 
while the family of unistochastic matrices B = B(a) forms the same hypocycloid in the two-dimensional cross-section 
of the four-dimensional body of TV = 3 bistochastic matrices determined by P,P 2 and P 3 . 

APPENDIX D: INTERPOLATION BETWEEN CORNERS OF TWO HYPOCYCLOIDS 



In this appendix we provide a discuss a family of unistochastic matrices of size 4, the spectra of which are not 
contained in the sum of 2, 3 and 4-hypocycloids. To find an interpolation between the corners of 3 and 4-hypocycloids 
consider the orthogonal matrix Os^cp) = 04(^)^1234, as defined in section IV, 



10 

10 

cos ip shop 

— sin if cos ip 





- 


1 





- 












1 


















1 






. 1 








. 





10 

cos (p sin ip 

— sin tp cos ip 
10 



(Dl) 



For ip varying in [0, tt/2] this family interpolates between a four-elements permutation P1234 and a matrix, the absolute 
values of which represent a three-elements permutation Pi24,3- Thus the spectra of the corresponding bistochastic 
matrices give an interpolation between the third and the fourth roots of identity. As shown in Fig. 4 this interpolation 
is located outside the hypocycloids H 3 and H4, so the support T,^ is larger than their sum G4. Repeating the 
argument with the multiplicative interpolation between any such a matrix and the Fourier matrix F 1 -- 4 ^ we conclude 
that all points inside the set bounded by this interpolation belong to the support An analogous scheme allows 
us to find an (N — 1) < — ► N interpolation. For example, the family of orthogonal matrices 0^5 ((f)) — 05(^)^12345 of 
size 5 gives an interpolation between l 1 / 5 and i = l 1 / 4 . 

To find the missing TV = 4 interpolation for the negative real part of the eigenvalues consider a permutation of the 
orthogonal matrix which contains the block O3 responsible for the 3-hypocycloid, 



O 



3,2 



" 1 








" 




" 


1 





0" 




" 


1 





" 





a 


b 


c 




1 













a 





c 


6 





c 


a 


b 













1 




c 





b 


a 


. 


b 


c 






. 





1 


. 




. b 





a 


c _ 



(D2) 



where its elements a, 6, and c are function of the angle ip as given in (4.11). Then the spectra of the corresponding 
orthostochastic matrices, obtained by squaring the elements of the above orthogonal matrices, 
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#3,4 (¥>) 



10 

cos 2 <p sin 2 if 

sin 2 ip cos 2 ip 
10 



and B 3i2 (ip) 



10 

a 2 c 2 b 2 



c 2 b 2 
b 2 a 2 



(D3) 



provide the required interpolations located outside hypocycloids H3 and H4 (see Fig. 4b and 4d). For larger dimen- 
sionality an analogous construction has to be performed to get the interpolations between the neighbouring roots of 
identity. 



APPENDIX E: AVERAGE TRACES OF UNISTOCHASTIC MATRICES 



To evaluate the average traces we rely on the results of Mello [M , who computed various averages, (.) , over the Haar 
measure on the unitary group U(N). In particular, he found average values of the following quantities constructed of 
elements U a b of a unitary matrix of size iV 



QlT-'-ZZ ((^ia • • • u b m pJiu aiai ■ ■ ■ u akak y). (ei) 

Mean trace of a squared unistochastic matrix, defined by By = |f/jj| 2 , reads 



N N 

t N , 2 := (TtB 2 ) use = ( £ (U lk U; k )(U kl U* kl )) u{N) = £ Qjjg = NQ^\ + N(N - 1)Q%£, (E2) 

i,k—l i,k=X 

since the symmetry of the problem allowed us to group together the terms according to the number of different indices. 
Using the results of Mello Q\\'^\ = 2/(N(N+l)) and Q^ 2 2 \ = l/(N(N + l)) we get t N . 2 = (N + 2)/(N + 1). 
The mean trace of B 3 reads 



tNy 



(TtB 3 



'USE 



NQi 



11,11,11 
1,11,11 



3N(N 



WniHl+N(N-l)(N-2)Q 1 1 



12,23,31 
2,23,31' 



(E3) 



The data in Mello's paper allow us to find Qxi'n'n = 6/[JV(JV 



1)(7V + 2)],Q™ 



1/[{N + 2){N 2 - 1)], and 

Qi^S^i = ~ 2)/[(N(N 2 - l)(N 2 - 4)], where the last term (wit h th ree different indices) is present only for 
N > 3. Substituting these averages into (E3) we arrive with the result ( 4.21 ). 
In the general case of arbitrary n we may write a formula 



tN': 



(TrB n ; 



USE 



i2,i2»3,- 



-Hn 
-lin 



(E4) 



which is explicit, but not easy to simplify. An analogous computation for the ensemble of sy mm etric unistochastic 
matrices may be based on results of Brouwer and Beenakker |57|, who computed the averages (El) for COE. 
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